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I.  INTRODUCTION 


In  general,  when  a  material  layer  is  accelerated 
by  a  light  pusher  such  as  provided  by  magnetic  pressure^ 
the  flow  is  unstable  to  the  Rayleigh-Taylor  instability 
at  the  interfaces.  A  perturbed  uniform  layer  rapidly 
forms  a  spike-bubble  structure  and  may  even  rupture. 

It  has  been  shown  previously  (Verdon  et  al.  1982)  that 
the  basic  dynamics  of  this  instability  process  in 
real,  incompressible  fluids  can  be  understood  on  the 
basis  of  the  irrotational  flow  of  layers  of  ideal, 
incompressible  fluids. 


In  Section  III,  we  demonstrate  the  equivalence  of  a  time- 
dependent  pressure  drive  and  a  gravity  drive  for  planar 
geometries.  We  also  discuss  the  dynamics  of  plane  thin 
shells  and  thin  shells  in  cylindrical  geometries.  In 

Section  IV,  the  mathematics  of  axisymmetric  flow  are 
discussed. 


II.  MATHEMATICAL  FORMULATION 

The  irrotational  flow  of  an  incompressible, 
inviscid  fluid  with  constant  density  is  completely 
specified  by  knowledge  of  the  velocity  potential  (ti 
which  must  satisfy  Laplace's  equation  wherever  the 
fluid  lies.  In  order  to  determine  ij)  ,  conditions  must 
be  given  at  the  surfaces  of  the  fluid.  Bernoulli's 
equation  may  be  used  as  an  evolution  equation  to  specify 
4>  at  the  surfaces.  Thus  the  mathematical  problem  is 
to  solve  Laplace's  equation  with  Dirichlet  conditions  on 
complicated  boundaries.  There  are  several  boundary 
integral  formulations  that  solve  such  problems;  in 
particular,  source  or  dipole  distributions  may  be 
used  which  lead  to  Rredholm  integral  equations  of  the 
first  or  second  kind  respectively.  The  mathematical 
properties  of  Fredholm  integral  equations  of  the  second 
kind  guarantee  efficient  numerical  techniques  for  their 
solution,  thus  providing  a  decided  advantage  to  the  use 
of  dipole  distributions.  Of  course,  one  cannot  solve 
all  potential  problems  via  dipole  distributions  alone; 
there  are  occasions  where  contributions  from  a  source 
or  sink  must  be  included,  but  this  step  is  relatively 
straightforward. 

Suppose  that  an  incompressible,  inviscid  fluid 
of  density  P  lies  between  two  surfaces  and  $2  • 
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Adjacent  to  the  fluid  layer,  there  lies  fluid  so  low 
in  density  that  the  region  it  occupies  may  be  considered 
effectively  a  vacuum,  yet  it  is  capable  of  supporting  an 
external  pressure.  Alternatively  the  region  is  a 
vacuum  but  there  are  external  surface  forces  present  such 
as  a  J  B  force  where  j  is  restricted  to  lie  in  the 
surface.  In  either  case,  assume  that  the  external  pressure 
is  constant,  but  with  different  values  Pj^  and  P2  in 
each  external  region  adjacent  to  Xj^  and  X2  »  respectively  < 
see  Figure  1  for  a  schematic  of  the  flow  geometry.  Clearly 
the  pressure  difference  will  accelerate  the  fluid  layer, 
inducing  a  mean  flow. 

Thus  the  velocity  potential  may  be  written  as 


(1) 


where  4>_  accounts  for  the  mean  (external)  flow  and 

c 

♦jj  may  be  expressed  in  terms  of  dipole  distributions 
along  each  surface  s^ 


^^(P) 


I 

I  j  u.  (q)  n(q)  •  V  G(p,q)dq 

=  1  c  *  ^ 


(2) 


where  p,q  are  field  points,  G  is  the  free  space 

A 

Green's  function,  n  is  the  normal  as  shown  in  Figure 

1,  and  is  the  gradient  operator  with  respect  to 

q.  The  choice  for  depends  on  the  geometrical 

£• 
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configuration  of  the  layer;  specific  cases  will  be 
discussed  later. 


The  velocity  u  of  the  surfaces  may  be  calculated 
from  u  =  7()i  .  In  particular,  the  velocity  potential 
takes  on  the  values  at  the  surfaces  s^ 

respectively,  where 


-  j  U  2^?)  (PCS2)  C4) 


Tangential  derivatives  of  together  with  the  external 

flow  give  the  tangential  velocities  of  the  surface, 
but  the  normal  component  obtained  by  differentiating 
(2)  involves  an  awkward  integral  to  be  performed  and 
it  is  preferable  to  proceed  via  a  different  approach. 

Introduce  the  vector  potential  X 

u  =  7  *  X,  with  7  •  X  =  0 
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(5) 


In  terms  of  dipole  distributions 


A(p)  =  I  /  U.{q)n(q)’<  V  G(p,q)dq 

£ _ *1  ^  i  y 


(6) 


1=1  s^ 


and  the  normal  velocity  follows  from  n  •  (V  A)  which 
involves  only  tangential  derivatives  of  A.  Thus,  once 
are  known,  the  surface  velocity  u^  may  be  computed 
as  described  above  and  the  surface  locations  updated 
in  time. 

So  far  the  kinematics  have  been  satisfied;  the 
dynamic  considerations  provide  an  evolution  equation 
for  y^.  The  starting  point  is  Bernoulli's  equation 
evaluated  at  the  surfaces. 


3(j, 


Di 


34) 


at 


E 

at 


I  (u^)^  Pi  “  (i=1.2)  (7) 


The  partial  time  derivative  is  Lagrangian  in  that  it 
represents  the  change  in  the  velocity  potential  following 
the  surface  motion.  The  substitution  of  (3)  and  (4) 
into  (7)  leads  directly  to  Fredholm  integral  equations 

ay^ 

of  the  second  kind  for  of  the  form; 


ay. 


3y. 


^(p)  +  2  f  — (q)n  (q)  .  V  G  (p,q)dq 


s . 
1 


ay. 


+  2  /  (q)n{q)*VgG(p,q)dq  =  )  (pcs,),  (8) 

S  M  “X 
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3y,  3w, 

(P)  +  2  /  (q)n(q) ‘V  G(p,q)dq 


+  2  ^  -j|  (q)n(q)-7qG(p,q)dq  =  ^2  1  ^2  '  ^2 '‘^e’ 


OMi 

where  contains  all  terms  without  • 

The  homogeneous  parts  of  equationa  (8)  and  (9) 
have  nontrivial  solutions  which  reflect  the  fact  the 
velocity  potential  is  determined  only  up  to  a  constant. 
According  to  the  Fredholm  alternative,  equations  (8) 
and  (9)  have  solutions  only  if  obey  a  certain 

condition.  Let  be  the  solutions  to  the  adjoint 
equations , 


Tj^(p)  -  2  f  Tj^(q)n(p)*V  G(p,q)dq 


2/  T2(q)n(p).VpG(p,q)dq  =  0,  (peg) 


T2(P)  +2  /  Tj^(q)n(p) -V  G{p,q)dq 

S  -  P 


2  f  2  ^  ^  =  0,  (P^  S- 


then  the  condition  for  a  solution  to  (8)  and  (9)  is 
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This  condition  actually  provides  a  relationship 
between  p^,  p^  and  the  external  flow  4)^.  In  the 
following  sections,  we  will  consider  various  geometries 
configurations  of  accelerating  thin  fluid  layers,  in 
which  case  the  form  of  4)_  will  be  known  since  it  mus 
describe  the  potential  flow  that  accelerates  the  fluid 
layers.  The  strength  of  the  acceleration  may  be 
directly  related  to  the  pressure  difference  p^  -  P2 
across  the  fluid  layer  by  means  of  (12). 


vr 


III.  PLANAR  LAYERS  AND  CYLINDRICAL  FLUTE  MODES 

For  planar,  two-dimensional  fluid  flow  (see 
Figure  la)  the  external  velocity  potential  is  obviously 


=  iv(t)y 


(13) 


where  v(t)  is  the  external,  uniform  velocity  of  the 
fluid  layer.  It  is  convenient  to  introduce  complex 
notation  and  to  let  2  =  x  +  iy  describe  a  field  point. 
The  surface  locations  may  then  be  parametrized  as 
Zj  (:t)  =  Xj  (a)  +  iy^  (a)  ,  j  =  1,2.  The  vector  potential, 

which  has  only  one  component  the  streamf unction  'i'  , 
may  be  combined  with  the  velocity  potential  into 
the  complex  potential  '?>  =  4>+  i'J'  .  In  particular 


=  -iv(t) 2  (14) 

E 

The  motion  of  the  fluid  surfaces  is  described 
by 


iv  (t) 


(15) 


where  the  star  superscript  implies  complex  conjugation. 
The  derivation  of  (15)  and  subsequent  equations  follow 
closely  the  approach  adopted  in  Baker  et  al  (1982).  The 
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3u  . 


Fredholm  integral  equations  for 


are 


_(a)  -  Re  f  -^(cf)  z-(cc)-z^(aM 


^la  )d0‘ 


i/- 


3u2{Q‘')  22^(ot')da' 


3t 


2i  (a)-z^  (a 


n) 


'  0l  •  2  ||  yi  (16) 


^‘'2,.  ....  I  1  ,*“l,  .. 


U<''>  *  hr 


Z2 (2)-z^  (a  ■  ) 


1  3Up  z-  (a*)da' 


-  12  -  2  (17) 


where  the  subscript  a  implies  differentiation  with 
respect  to  a  and 


=  Re 


i  f 

TT  1  ’ 


(a' ) 


z^  (a) -2^  (a' ) 


q  .  (a  '  ) 
^la  ' 


°1  "‘^1  ) 
z^  (a)  -2^  (a  '  ) 


U  2  ( '  ) 
(a)-Z2 (a' ) 


q^^  (°‘)-q2  ) 

z^  (a)-22  (a-  ) 


da' 


-  + 


2Pi 


(18) 
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92  *  “Re 


(a  • ) 


92 (a)“qi  («' ) 

(a)-Zj^  (a* )  *la^“  ^ 


U2(ci') 

Zj  (a)-22  (ex'  ) 


q,  (Qi)-q,  (a’ ) 

zJlaPI^TcPT 


+  9292  "  2P2 


(19) 


In  order  for  (17)  and  (18)  to  have  a  solution,  the 
condition  (12)  must  be  satisfied;  thus 


da  +  ! 


f  dot  -  /  da 

(20) 


For  given 


and  p 


2' 


(20)  determines 


dv 

3t- 


The  form 


of  the  equations  (16)  and  (17)  is  identical  to  those 
for  a  fluid  layer  falling  under  the  influence  of  gravity 


with  g  =  -  ^ 


dv 

dt 


The  above  analysis  demonstrates  the  equivalence  of 
pressure  drive  (normal  to  the  surface)  in  planar  geometries 
and  a  time-dependent  gravity  drive  (in  a  fixed  direction) . 
Eq.  (20)  guarantees  conservation  of  momentum  in  these 
systems  and  so  determines  the  mean  acceleration  of  the  layer. 

The  numerical  solution  to  (15),  (16),  and  (17) 

follows  a  standard  approach  adopted  by  Balter  et  al  (1982). 

We  present  here  the  results  for  ^  ~  The  initial 

dt 

profile  has  the  form 
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a  +  i 


22  =  a  +  i  0.06  ir  cos  a. 

Figure  2  shows  the  evolution  in  time  of  the  surface 
locations.  A  clearly  observed  spike  develops  and  the 
layer  thickness  of  the  bubble  thins  dramatically. 

Figure  3  shows  the  layer  thickness  at  the  bubble  center 
as  a  function  of  time.  Unfortunately , the  code  presently 
is  not  accurate- beyond  a  time  of  approximately  4.  Further 
calculations  in  time  require  a  modification  of  the  numerical 
quadrature  for  the  integrals  in  fl6)  ,  (17) ,  (18) ,  and  (19). 

Next,  we  consider  cylindrical  geometry  as  shown 
in  Figure  lb.  The  external  velocity  ootential  now 
has  the  form 

4i_  =  -  A(t)  log(z)  (21) 

£ 

An  equation  similar  to  (20)  is  easily  obtained,  but 
there  is  an  important  difference  between  the  previous 
case  and  the  present  one;  the  external  flow  no  longer 
has  the  appearance  of  that  induced  by  a  gravity  field. 
Numerical  calculations  give  the  results  shown  in  Figure  4. 
The  spike  and  bubble  structure  are  different  from  those 
observed  in  Figure  2. 

Finally,  in  the  next  section,  we  present  the 
mathematical  formulation  for  axisymmetric  flow,  which 
is  somewhat  different  from  that  considered  in  this  section. 
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IV.  AXISYMMETRIC  FLOW  AND  SAUSAGE  MODES 


The  generalized  vortex  method  of  Baker,  Meiron  and 

Orszag  (1980,1982)  is  a  boundary  integral  technique  applied 
to  free  surface  flow  in  incompressible,  inviscid  fluid 
which  contain  regions  of  constant,  but  different  densities. 
When  the  free  surface  lies  betv>re>^n  fluid  and  vacuum,  a 
simpler,  but  equivalent  metho!  may  be  used.  This  method 
is  described  in  some  generality  before  being  applied  to 
axisymmetric  flow. 

The  free  surface  is  represented  as  a  dipole  layer  of 
strength  u.  The  velocity  potential  <p  may  then  be 
written  as 

<P(P)  =  J  u(q) 

where  G  is  the  free  space  Green's  function  for  Laplace's 

equation,  p  and  q  are  field  points  and  n  is  an 
outward  normal.  For  convenience  we  assume  that  the  free 
surface  is  smooth  (has  a  continuous  normal)  and  closed; 
an  open  surface  may  be  considered  as  a  particular  limit 
of  a  closed  surface.  The  normal  derivative  of  G  is 
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taken  with  respect  to  q,  keeping  p  fixed.  As  p 
approaches  the  free  surface  along  the  normal,  the  potential 
takes  on  different  values  on  either  side 


<^j(P)  =  f  If  +  |(P)  (23) 

as  p  tends  to  the  surface  from  the  inside  and 

♦q(P)  =  f  p(q)  |f(P»q)dq  -  ^(p)  (24) 

as  p  tends  to  the  surface  from  the  outside.  Clearly, 
y  (P)  =  (pj  (P)  -<|>q(P)  . 

If  y (P)  is  known,  then  the  potential  can  be  evaluated 
from  (22)  and  the  fluid  velocity  follows  from  u  =  V(j,  . 
Alternatively  one  may  obtain  the  fluid  velocity  from  the 
velocity  vector  potential  a  ;  u  =  v  x  A.  For  a  dipole 
layer 

A{p)  =/  yfq)  nx  vG(P»q)  dq  ^25) 

where  the  gradient  is  taken  with  respect  to  q  and  h 
is  the  unit  outward  normal.  As  p  tends  to  the  surface, 
X(p)  is  continuous. 

Normally  in  free  surface  flow  problems,  the  velocity 
is  required  only  at  the  surface  in  order  to  update  its 
location.  This  is  easily  accomplished  if  g  is  known. 


1 

1 

In  fact,  the  tangential  velocity  components  can  be  computed 
from  V<(i  and  the  normal  component  from  V  A.  In  both 
cases,  only  tangential  derivatives  are  required  thus  avoiding 
the  more  complicated  and  more  difficult  evaluation  of  normal 
derivatives . 

The  basis  of  the  generalized  vortex  method  is  the  use 
of  an  evolution  equation  for  p  derived  from  Bernoulli's 
equation,  by  which  p  is  updated  and  the  free  surface 
velocity  calculated  as  described  above.  The  method  involves 
evaluating  the  time  derivatives  of  (23),  (24)  and  (25), 
which,  in  the  case  of  axisymmetric  flow,  is  tedious 
algebraically  and  presents  difficulties  numerically. 

However,  if  one  is  interested  in  the  flow  of 
a  free  surface  between  fluid  and  vacuum  one  may  use 
Bernoulli's  equation  directly  to  update  the  potential  at 
the  free  surface  as  it  moves.  Then  (23)  or  (24)  is  used 
as  a  Fredholm  integral  equation  of  the  second  kind  for  p  . 

Finally,  knowing  p,A  is  computed  from  (25)  and  the  free 
surface  velocity  is  computed  as  described  above.  This  is  the 
approach  adopted  for  axisymmetric  free  surface  flow. 

Introduce  a  cylindrical  coordinate  system  (R,9,z)  where 
R  is  the  radius,  6  the  aximuthal  angle  and  z  the  axial 
position.  The  free  surface  is  pareimetrized  by  R(e,t)  and 
z(e,t).  Let  1(e)  be  the  principal-value  integral  in 


1(e)  .  ^  # 


M(e»)de‘ 

[(z(e)-2(e‘)) 

Zg  (e ' ) [ (2 (e) -z (e ' ) ) (e) -R^ (e ' ) 1 -2R(e • ) Rg (e ' ) (z (e) -z (e • ) ) 

(z(e)-z(e'))^+(R(e)-R(e'))^ 

X  E(k)|  (26) 

where  K(k)  and  E(k)  are  the  complete  elliptic  integrals 
of  the  first  and  second  kinds  respectively,  the  subscript 
e  refers  to  differentiation  with  respect  to  e  (t  fixed) 
and 


^+(R(e)+R(e'))^l^''^ 


2g(e')K(k) 


j^2  ^  _ 4R(e)R(eM _ 

(z  (e)-2 (e') )^+(R{e)+R(e' ) ) 


(27) 


The  only  non-vanishing  component  of  X  is  the  azimuthal 
component,  expressed  as  i|</R  ,  where  plays  the  role 
of  the  streamfuncti on  in  axisymmetric  flow.  At  the  surface, 


(j)  (e) 


u (e* )de' 
"2^ 


^^[(z(e)-z(e’))^+(R(e)+R(e'))^]^^^ 


{  ^*(®>"2(e’))-R(e')R_(e')]K(k) 


-(Zg{e’ )  (z  (e)-z  (e* )  )-R(e’ )  Rg(e* )  ] 


(z(e)-z(e')  )^+R^(e)+R^(e')  ]+2R^(e)R(e')Rg(e') 

- 2 - 2 - 

(z(e)-z(e*) )^+  (R(e)-R(e')) 
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(28) 


The  radial  and  axial  velocity  components,  u  and  w 


respectively,  are  given  by 

n  -  ii  -  1  ii 

3R  “  R  9z 


(29a) 


w  =  11  =  -  i  M 
3z  R  3R 


(29b) 


TO  update  the  free  surface,  the  velocity  components  at  the 
surface  must  be  known.  They  can  be  expressed  in  terms  of 
the  derivatives  of  (P  and  ip  along  the  surface.  Using 
(29)  one  finds 


R  11  +  z  11  = 

e  3R  e  3z 


R^u  +  z  w 
e  e 


(30a) 


li  +  2 

3R  e 


3z 


R(z  u 
e 


R  w) 
e 


(30b) 


where  (|)  will  either  be  or  depending  on  which 
side  of  the  interface  the  fluid  lies.  Inverting  the 
above  equations,  one  obtains 


3R  _ 
3t 


u 


(31a) 


32 

3t 


w 


<Ve  - 


Ve/«)^e 


(31b) 
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as  the  evolution  equations  for  the  surface's  location, 
where 


s 


2 

e 


2  2 

K  + 

e  e 


(32) 


Finally,  Bernoulli's  equation  is  used  as  an  evolution 
equation  for  41  .  As  a  specific  example,  consider  a 
spherical  bubble  that  is  compressed  by  an  external  pressure 
gradient  in  the  fluid.  Write  the  potential  as 


A(t) 


+  4) 


(33) 


where  A  is  the  strength  at  a  sink  (or  source,  if  A<  0) 
placed  at  center  of  the  bubble.  After  substituting  (33) 
into  Bernoulli's  equation,  one  finds  at  the  surface 


34i  _ 
at 


dA 

dt 


4 


— 5 — TT 


+  Pi 


Po 


(34) 


where  the  time  derivative  is  taken  following  the  motion 

described  by  (31)  the  fluid  is  assumed  stationary  far 

from  the  bubble  with  pressure  ,  the  density  is  1  and 

P^  is  the  pressure  inside  the  bubble.  If  the  bubble  is 

2  2  2 

a  perfect  sphere,  one  may  set  ^  *  0  and  p  ■  R  +  2  so 
that  (31)  and  (34)  become 


(35) 


1  dA  ,  A _  _ 

p  0  I 


dt 


"  ^ 


(36) 


We  may  subtract  (35)  from  (34)  to  obtain 


_  i(u2+w2)  +  ^ 

3t  2^“  '  dt 


(37) 


Also, 


3R  =  n  -  AR 

3^ 


(38a) 


AZ 


(38b) 


where  u  and  w  are  computed  from  the  dipole  representation 
for  as  described.  The  above  equations,  (35),  (36),  (37) 
emd  (38)  constitute  a  set  of  evolution  equations  for  the 
motion  of  the  bubble  surface. 

The  method  of  solution  is  as  follows.  Suppose 
A,  p,  R,  z  and  ^  are  all  known  at  some  given  time.  The 

time  derivatives  of  A  and  p  are  known  from  (35)  and 
(36).  We  then  solve  for  u  by  iteration  from  (24) 
with  ♦q  ■ 
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»  -  2<|»(e)  +  2l^(e) 


(39) 


where  represents  1(e)  with  pCe')  replaced  by 

Ujj(e').  The  iteration  process  is  known  to  be  globally 
convergent.  Moreover  standard  extrapolation  techniques 
may  be  used  to  ensure  that  only  a  few  (2  ~3)  iterations 
are  required  for  accurate  results.  Knowing  p  (e)  ,  ^j(e) 

is  computed  from  (28)  and  then  u  and  w  are  determined 
from  (31).  Finally,  the  time  derivatives  of  R  and 
z  are  evaluated  by  (37)  and  (38)  and  the  surface  may 

be  updated. 

To  execute  this  method  of  solution  ntunerically  requires 
some  care.  Points  evenly  distributed  in  e  are  introduced 
as  a  discrete  representation  for  the  surface.  All 
derivatives  in  e  are  approximated  by  spline  derivatives. 

The  difficulties  lie  in  the  numerical  evaluation  of  the 
integrals  (26)  and  (28).  The  complete  elliptic  integral 
of  the  first  kind,  K(k),  has  a  logarithmic  singularity 
as  k  1  which  occurs  when  e'-*-e.  The  complete 

elliptic  integral  of  the  second  kind  behaves  as 
(e'-e) log] e'-e| 

as  €'-►  e.  A  further  complication  is  the  fact  that,  for 
e  close  to  the  axis,  k  varies  from  0  at  the  axis  to 
1  at  e'  ■  e.  This  rapid  variation  in  k  must  be  accounted 
for  in  any  accurate  numerical  quadrature. 

The  integrals  may  be  regularized  by  using  the  fact 
that,  when  v  ■  1,  1(e)  =  y  and  i);(e)  ®  0.  In  addition, 
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a  particular  quadrature  has  been  devised  to  account  for 
the  presence  of  terms  containing  logarithms.  Suppose  that 
the  parametrization  of  the  free  surface  is  such  that 
0  ^  e  and  R(0)  *  RCif)  »  o,  that  is  e  *  0,  if 

correspond  to  the  poles  of  the  bubble.  K(k)  and  E(k) 
may  be  approximated  by 

K(k)  »  Aj^(m)  -  Bj^(m)  log  |  m|  (40a) 

E(k)  =  Aj.(m)  -  Bg(m)  log  |  m]  (40b) 

(m)  ,  A^(m),  Bj^(m)  and  Bp(m) 

are  polynomials  in  m  (see  Abramowitz  and  Stegun,  1964). 
Thus  one  nay  split  the  integrals  into  parts  that  contain 
logarithms  or  not.  Define 


where  m  =  1  -  k  and 


F(A,B)  = 


u (e’ )  -  u (e) 


((z(e)-z(e’))^+(R(e)+R(e'))^)^^^ 


jzg(e'). 


(e-)  [(z(e)-z(e'))^-*-R^(e)V(e’)]-2R(e’)Re(e')  (2(e)-z(e«))  g  | 

(2(e)-z(e'))^+(R(e)”B(®*)) 

(41) 


II  «  ^  I 


(42) 
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I2  =  ^  /  I ‘’I  ^ 

^3  "  W  I  1"^!  (^) 


Then  I  =  -  I2  ~  I3.  Here  q  models  the  behaviour  of 

m;  m  =  1  at  e*  =  0,ti  and  m  =  0  at  e'  =  e,  so 
the  variation  in  m  is  rapid  when  e  is  near  0  or  tt  , 
A  good  choice  for  q  is 


q 


e'  (Tr-e’)+o  (e'-e)^ 


(45) 


where  0  =  lA  ^  +  tt  [4e  (tt -e)  ]  .  Since  F(A^,Ag)  and 

F(B„,B„)  are  analytic,  I,  and  I,  may  be  evaluated  by 

2 

the  trapezoidal  rule  with  0 (h  )  accuracy , where  h  is 
the  spacing  in  e  of  the  points.  For  a  special 

quadrature  is  used.  P(B„,B„)  is  approximated  as  a 
piecewise  continuous  linear  function  and  is  then 

integrated  analytically;  this  leads  to 

I3  =  I  ‘"j  (46) 

where  F.  is  the  value  of  F(B„,B_)  at  e'..  After 
3  R  fc  j 

some  algebra,  one  finds 
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(47) 

where 


Aj  =  j(e^-ej\og  {o  (ej-e)^} 


2  (o  -1)  ^  (e'-e)  ^-r^+2s  (o-l) 
4  (a-1)^ 


log|(o-l) (e!-e) ^+r (e!-e) +s  } 


^  r  [4s(o-l)-r^l 
2(0-1)  ^  (0-1)^ 


1/2 


tan 


1 


2 (o-l) (el-e )+r 
(4s(a-l)-r^l^/^ 


Bj  -  (e^-e)log  {o  (ej-e)^  ) 


(48) 


2  (o -1)  (e!-e)+r  2 

- — — rJ - log  {  (o-l)  (e!-e)  +r(e!-e)+s  } 

2  ( 0-i )  J  J 


(4s(a-l)-r“l^/^  .  -1 

-  tan 

(0-1) 


j  ? ( 0-1) (e^-e) +r 


[4s(a-l)-r‘l 


27171 


(49) 


and  r  =TT-  2e,  s  =  e(Tr-e). 

The  total  quadrature  error  for  I  is  of  order  O(h^) 
near  the  axis  and  is  more  accurate  away  from  the  axis.  A 
quadrature  for  i(i(e)  is  similarly  constructed;  define 
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Q  _  [p  (e*)  -  M  (e)  ] _ 

t  (2(e)-2(e')  )^+(R(e)+R(e'))^l 


Iz-(e')  (2(e)-2(e'))-R(e')R^(e')]  (50) 


tu  (e* )  -  u  (e) 


[ (2 (e)-z(e*) / ^+(R(2)+R(e' ) 


jlZg(e')  (2(e)-z(e’))-RteMReteM] 


((z(e)-z(e'))^+R^(e)+R^(e')]  +  2R^  (e)  R(e ' )  Rg  (e '  )j  (53^) 

*1  = 


.(z  (e)-z(e’)  )^+(R(e)-R(e 


- 2“  ”2""“^^ - 2I  I 

•))2  sfce-eM^J  ^ 


_  1  ",  H 

'2  ■  ^  (,  s2(e--e)  * 


*^3  ■  5F 


H  B  (m) 

- -= - 5-  {log  lin|-  log  jqj  }  ]  de' 

(2  (e)-2(e'))‘'+CRCe)-R(e')) 


i 


r,  H  B  (m) 

I  - _ 

0  (2(e)-z(e'))^+(R(e)-R(e'))^ 


log  I q I  de ' 


(55) 


Then  'I'  *  ~  ii'2  "  ’^'3  *  'l’4‘  integrals  and 

are  evaluated  by  the  trapezoidal  rule,  and  is 

evaluated  by  the  special  quadrature  (46) .  However 
has  a  pole  singularity,  since  H  - (e'~e)  as  e' 
approaches  e.  We  devise  a  quadrature  by  approximating 
H  =  HBg(m)/tSg(e'-e) ]  as  a  piecewise  continuous  linear 
function  and  integrating  analytically;  thus 

ij;-  =  ^  w.  H.  (56) 

j  ^  ^ 


where 


“3  '  K 

+  (e'_^-e)log|e!_^-e|-2(e*-e)log|e--e|)  (57) 

Starting  with  the  initial  conditions  A  =  0,  ^  =  0, 

R=  (1  ■*"  ccos(2e)  )sin(e)  ,  z  «  -  (1  +  e  cos  (2e)  )cos  (e)  , 
p  *  0.9703797,  the  surface  is  discretised  with  33  points 
and  updated  in  time  with  a  fourth  order  Adams-Bashforth 
method  with  a  timestep  of  0,01.  The  pressure  difference 
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r 


Pq  -  Pj  =  1  and  e  =  0.1.  Figure  5  shows  the 
profile  of  the  bubble  at  various  times.  Clearly  a  Rayleigh- 
Taylor-type  instability  develops  at  the  waist  of  the  bubble 
and  will  probably  lead  to  bubble  separation  into  two 
smaller  bubbles. 
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AA  =  a05 


(d)  Interface  geometry  at  t  =  1.396. 


CYCLE  =  570  TIME  =  3.979 


(c)  Interface  at  t  =  3.979. 
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Figure  3.  A  plot  of  shell  thickness  versus  time  for 
the  flute  mode  instability  problem. 
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(b)  Outside  of  a  shell  only 
one- third  as  thick  (kAx 


(a)  Outside  of  a  relatively 
thick  shell  (kAx  =  3). 

Figure  4.  Results  of  a  calculation  with  a  sinall 
perturbation  (0.57o). 


